/**********************************************************************//**
 @file primitives.hpp 
 @brief Contains some primitives datatypes and methods for computational geometry.
 @author Ricardo C. Marques <rmarques@inf.puc-rio.br>
 @version 0.0.1
 @date 05/11/2012
**************************************************************************/
// Standard Libraries
#include <stdlib.h>
#include <math.h>

// Module Own Header File
#include "primitives.hpp"

// Auxiliary Header Files
extern "C" 
{
    #include "predicates.h"
};

double orient2d(Point2D A, Point2D B, Point2D C)
{
    return orient2d((double*)&A,(double*)&B,(double*)&C);
}
double orient2d(Point3D A, Point3D B, Point3D C)
{
    return orient2d((double*)&A,(double*)&B,(double*)&C);
}
double incircle(Point2D A, Point2D B, Point2D C, Point2D D)
{
    return incircle((double*)&A,(double*)&B,(double*)&C,(double*)&D);
}
double incircle(Point3D A, Point3D B, Point3D C, Point3D D)
{
    return incircle((double*)&A,(double*)&B,(double*)&C,(double*)&D);
}
bool OrientedPositively(Point2D A, Point2D B, Point2D C)
{
    return orient2d(A,B,C)>0.;
}
bool OrientedPositively(Point3D A, Point3D B, Point3D C)
{
    return orient2d(A,B,C)>0.;
}
static bool isCoincident(Point2D A, Point2D B)
{
    return  (A.x==B.x)&&(A.y==B.y);
}
static void Add(Point2D A, Point2D B, double C[2])
{
    double ZeroZero   [2]={ 0.,0.};
    double MinusOneOne[2]={-1.,1.};

    double AxBx[2]={A.x,B.x};
    double AyBy[2]={A.y,B.y};
    C[0] = orient2d(AxBx,MinusOneOne,ZeroZero); // Cx=Ax+Bx
    C[1] = orient2d(AyBy,MinusOneOne,ZeroZero); // Cy=Ay+By
}
static void Lerp(Point2D A, Point2D B, double t, Point2D* P)
{
    double ZeroA[2]={ 0., A.x};
    double tZero[2]={ t ,  0.};
    double OneB [2]={ 1., B.x};
    P->x=orient2d(tZero, OneB, ZeroA);
    ZeroA[1]=A.y; // ZeroA[2]={ 0., A.y};
    OneB [1]=B.y; // OneB [2]={ 1., B.y};
    P->y=orient2d(tZero, OneB, ZeroA);
}
static void Subtract(Point2D A, Point2D B, double C[2])
{
    double ZeroZero[2]={0.,0.};
    double   OneOne[2]={1.,1.};

    double AxBx[2]={A.x,B.x};
    double AyBy[2]={A.y,B.y};
    C[0] = orient2d(AxBx,OneOne,ZeroZero); // Cx=Ax-Bx
    C[1] = orient2d(AyBy,OneOne,ZeroZero); // Cy=Ay-By
}
static double Subtract(double a, double b)
{
    double ZeroZero[2]={0.,0.};
    double   OneOne[2]={1.,1.};
    double       ab[2]={a ,b };
    return orient2d(ab,OneOne,ZeroZero);
}
Point2D GetCircumcenter(Point2D A, Point2D B, Point2D C )
{
    Point2D O={0.,0.}; // circumcenter

    double AB[2];// = B-A
    Subtract(B,A,AB); 

    double AC[2];// = C-A
    Subtract(C,A,AC);

    double ZeroZero[2]={0.,0.};
    double delta=orient2d(AB,AC,ZeroZero);// = (Bx-Ax)*(Cy-Ay)-(By-Ay)*(Cx-Ax)
    if( delta==0. )
    {
        // Triangle is collapsed (with area=0)
        if( isCoincident(A,B) )
        {
            if( isCoincident(B,C) )// All 3 points are coincident
                return A;
            else
            {
                // A/Bo-------------oC
                Lerp(A,C,0.5, &O);
                return  O;
            }
        }
        else if( isCoincident(B,C) ) // C cannot be coincident to A
        {
            // Ao-------------oB/C
            Lerp(A,C,0.5, &O);
            return  O;
        }
        else if( isCoincident(A,C) )
        {
            // Bo-------------oA/C
            Lerp(A,B,0.5, &O);
            return  O;
        }
#if _DEBUG
        //else // Impossible case
            return  O;
#endif
    }

    double  AxAy[2]={ A.x, A.y};
    double mAyAx[2]={-A.y, A.x};
    double SqrLenA = orient2d(AxAy, mAyAx, ZeroZero);// = Ax*Ax + Ay*Ay

    double  BxBy[2]={ B.x, B.y};
    double mByBx[2]={-B.y, B.x};
    double SqrLenB = orient2d(BxBy, mByBx, ZeroZero);// = Bx*Bx + By*By

    double  CxCy[2]={ C.x, C.y};
    double mCyCx[2]={-C.y, C.x};
    double SqrLenC = orient2d(CxCy, mCyCx, ZeroZero);// = Cx*Cx + Cy*Cy

    double DiffSqrLenAB = Subtract(SqrLenB, SqrLenA);
    double DiffSqrLenAC = Subtract(SqrLenC, SqrLenA);

    double den = 2.*delta;
    double ABxACx[2]={AB[0], AC[0]};    
    double AByACy[2]={AB[1], AC[1]};
    double dABdAC[2]={DiffSqrLenAB, DiffSqrLenAC};
    O.x =  orient2d(AByACy, dABdAC, ZeroZero)/den; //num = -AC[1]*DiffSqrLenAB +  AB[1]*DiffSqrLenAC
    O.y = -orient2d(ABxACx, dABdAC, ZeroZero)/den; //num =  AC[0]*DiffSqrLenAB + -AB[0]*DiffSqrLenAC
    return O;
}
Point2D GetCircumcenter(Point3D _A, Point3D _B, Point3D _C )
{
    Point2D A = {_A.x, _A.y};
    Point2D B = {_B.x, _B.y};
    Point2D C = {_C.x, _C.y};
    return GetCircumcenter(A,B,C);
}